# Replication file for "Value Shift: Immigration Attitudes and the Sociocultural Divide" by Caroline Marie Lancaster. Please contact cml93@live.unc.edu with any questions. 

# This file compiles the LISS analysis dataset and creates figures 1 and 7.

# 3 August 2020

library(readstata13)
w1 <- read.dta13("cv08a_1.1p_EN.dta")
w2 <- read.dta13("cv09b_2.1p_EN.dta")
w3 <- read.dta13("cv10c_EN_1.0p.dta")
w4 <- read.dta13("cv11d_EN_1.0p.dta")
w5 <- read.dta13("cv12e_EN_1.0p.dta")
w6 <- read.dta13("cv13f_EN_1.0p.dta")
w7 <- read.dta13("cv14g_EN_1.0p.dta")
w8 <- read.dta13("cv16h_EN_1.0p.dta")
w9 <- read.dta13("cv17i_EN_1.0p.dta")
w10 <- read.dta13("cv18j_EN_1.0p.dta")
w11 <- read.dta13("cv19k_EN_1.0p.dta")
abortion <- read.dta13("Abortion.dta")
religion <- read.dta13("Religion.dta")
dem <- read.dta13('Demographics.dta')

#selecting variables from each Politics and Values dataset
w1 <- w1[c('nomem_encr', 'cv08a101', 'cv08a103','cv08a104', 'cv08a105', 'cv08a109', 'cv08a110', 'cv08a111', 'cv08a112', 'cv08a113', 'cv08a114', 'cv08a115', 'cv08a116', 'cv08a117', 'cv08a118', 'cv08a119', 'cv08a120', 'cv08a121', 'cv08a122', 'cv08a123', 'cv08a124', 'cv08a125', 'cv08a126', 'cv08a127', 'cv08a128', 'cv08a129', 'cv08a130', 'cv08a131', 'cv08a132', 'cv08a133', 'cv08a134', 'cv08a143', 'cv08a144', 'cv08a145', 'cv08a146', 'cv08a151', 'cv08a152', 'cv08a153', 'cv08a154')]

w2 <- w2[c('nomem_encr', 'cv09b101', 'cv09b103','cv09b104', 'cv09b105', 'cv09b109', 'cv09b110', 'cv09b111', 'cv09b112', 'cv09b113', 'cv09b114', 'cv09b115', 'cv09b116', 'cv09b117', 'cv09b118', 'cv09b119', 'cv09b120', 'cv09b121', 'cv09b122', 'cv09b123', 'cv09b124', 'cv09b125', 'cv09b126', 'cv09b127', 'cv09b128', 'cv09b129', 'cv09b130', 'cv09b131', 'cv09b132', 'cv09b133', 'cv09b134', 'cv09b143', 'cv09b144', 'cv09b145', 'cv09b146', 'cv09b151', 'cv09b152', 'cv09b153', 'cv09b154')]

w3 <- w3[c('nomem_encr', 'cv10c101', 'cv10c103','cv10c104', 'cv10c105', 'cv10c109', 'cv10c110', 'cv10c111', 'cv10c112', 'cv10c113', 'cv10c114', 'cv10c115', 'cv10c116', 'cv10c117', 'cv10c118', 'cv10c119', 'cv10c120', 'cv10c121', 'cv10c122', 'cv10c123', 'cv10c124', 'cv10c125', 'cv10c126', 'cv10c127', 'cv10c128', 'cv10c129', 'cv10c130', 'cv10c131', 'cv10c132', 'cv10c133', 'cv10c134', 'cv10c143', 'cv10c144', 'cv10c145', 'cv10c146', 'cv10c151', 'cv10c152', 'cv10c153', 'cv10c154')]

w4 <- w4[c('nomem_encr', 'cv11d101', 'cv11d103','cv11d104', 'cv11d105', 'cv11d109', 'cv11d110', 'cv11d111', 'cv11d112', 'cv11d113', 'cv11d114', 'cv11d115', 'cv11d116', 'cv11d117', 'cv11d118', 'cv11d119', 'cv11d120', 'cv11d121', 'cv11d122', 'cv11d123', 'cv11d124', 'cv11d125', 'cv11d126', 'cv11d127', 'cv11d128', 'cv11d129', 'cv11d130', 'cv11d131', 'cv11d132', 'cv11d133', 'cv11d134', 'cv11d143', 'cv11d144', 'cv11d145', 'cv11d146', 'cv11d151', 'cv11d152', 'cv11d153', 'cv11d154')]

w5 <- w5[c('nomem_encr', 'cv12e101', 'cv12e103','cv12e104', 'cv12e105', 'cv12e109', 'cv12e110', 'cv12e111', 'cv12e112', 'cv12e113', 'cv12e114', 'cv12e115', 'cv12e116', 'cv12e117', 'cv12e118', 'cv12e119', 'cv12e120', 'cv12e121', 'cv12e122', 'cv12e123', 'cv12e124', 'cv12e125', 'cv12e126', 'cv12e127', 'cv12e128', 'cv12e129', 'cv12e130', 'cv12e131', 'cv12e132', 'cv12e133', 'cv12e134', 'cv12e143', 'cv12e144', 'cv12e145', 'cv12e146', 'cv12e151', 'cv12e152', 'cv12e153', 'cv12e154')]

w6 <- w6[c('nomem_encr', 'cv13f101', 'cv13f103','cv13f104', 'cv13f105', 'cv13f109', 'cv13f110', 'cv13f111', 'cv13f112', 'cv13f113', 'cv13f114', 'cv13f115', 'cv13f116', 'cv13f117', 'cv13f118', 'cv13f119', 'cv13f120', 'cv13f121', 'cv13f122', 'cv13f123', 'cv13f124', 'cv13f125', 'cv13f126', 'cv13f127', 'cv13f128', 'cv13f129', 'cv13f130', 'cv13f131', 'cv13f132', 'cv13f133', 'cv13f134', 'cv13f143', 'cv13f144', 'cv13f145', 'cv13f146', 'cv13f151', 'cv13f152', 'cv13f153', 'cv13f154')]

w7 <- w7[c('nomem_encr', 'cv14g101', 'cv14g103','cv14g104', 'cv14g105', 'cv14g109', 'cv14g110', 'cv14g111', 'cv14g112', 'cv14g113', 'cv14g114', 'cv14g115', 'cv14g116', 'cv14g117', 'cv14g118', 'cv14g119', 'cv14g120', 'cv14g121', 'cv14g122', 'cv14g123', 'cv14g124', 'cv14g125', 'cv14g126', 'cv14g127', 'cv14g128', 'cv14g129', 'cv14g130', 'cv14g131', 'cv14g132', 'cv14g133', 'cv14g134', 'cv14g143', 'cv14g144', 'cv14g145', 'cv14g146', 'cv14g151', 'cv14g152', 'cv14g153', 'cv14g154')]

w8 <- w8[c('nomem_encr', 'cv16h101', 'cv16h103','cv16h104', 'cv16h105', 'cv16h109', 'cv16h110', 'cv16h111', 'cv16h112', 'cv16h113', 'cv16h114', 'cv16h115', 'cv16h116', 'cv16h117', 'cv16h118', 'cv16h119', 'cv16h120', 'cv16h121', 'cv16h122', 'cv16h123', 'cv16h124', 'cv16h125', 'cv16h126', 'cv16h127', 'cv16h128', 'cv16h129', 'cv16h130', 'cv16h131', 'cv16h132', 'cv16h133', 'cv16h134', 'cv16h143', 'cv16h144', 'cv16h145', 'cv16h146', 'cv16h151', 'cv16h152', 'cv16h153', 'cv16h154')]

w9 <- w9[c('nomem_encr', 'cv17i101', 'cv17i103','cv17i104', 'cv17i105', 'cv17i109', 'cv17i110', 'cv17i111', 'cv17i112', 'cv17i113', 'cv17i114', 'cv17i115', 'cv17i116', 'cv17i117', 'cv17i118', 'cv17i119', 'cv17i120', 'cv17i121', 'cv17i122', 'cv17i123', 'cv17i124', 'cv17i125', 'cv17i126', 'cv17i127', 'cv17i128', 'cv17i129', 'cv17i130', 'cv17i131', 'cv17i132', 'cv17i133', 'cv17i134', 'cv17i143', 'cv17i144', 'cv17i145', 'cv17i146', 'cv17i151', 'cv17i152', 'cv17i153', 'cv17i154')]

w10 <- w10[c('nomem_encr', 'cv18j101', 'cv18j103','cv18j104', 'cv18j105', 'cv18j109', 'cv18j110', 'cv18j111', 'cv18j112', 'cv18j113', 'cv18j114', 'cv18j115', 'cv18j116', 'cv18j117', 'cv18j118', 'cv18j119', 'cv18j120', 'cv18j121', 'cv18j122', 'cv18j123', 'cv18j124', 'cv18j125', 'cv18j126', 'cv18j127', 'cv18j128', 'cv18j129', 'cv18j130', 'cv18j131', 'cv18j132', 'cv18j133', 'cv18j134', 'cv18j143', 'cv18j144', 'cv18j145', 'cv18j146', 'cv18j151', 'cv18j152', 'cv18j153', 'cv18j154')]

w11 <- w11[c('nomem_encr', 'cv19k101', 'cv19k103','cv19k104', 'cv19k105', 'cv19k109', 'cv19k110', 'cv19k111', 'cv19k112', 'cv19k113', 'cv19k114', 'cv19k115', 'cv19k116', 'cv19k117', 'cv19k118', 'cv19k119', 'cv19k120', 'cv19k121', 'cv19k122', 'cv19k123', 'cv19k124', 'cv19k125', 'cv19k126', 'cv19k127', 'cv19k128', 'cv19k129', 'cv19k130', 'cv19k131', 'cv19k132', 'cv19k133', 'cv19k134', 'cv19k143', 'cv19k144', 'cv19k145', 'cv19k146', 'cv19k151', 'cv19k152', 'cv19k153', 'cv19k154')]

names(w1) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w1$wave <- 1

names(w2) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w2$wave <- 2

names(w3) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w3$wave <- 3

names(w4) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w4$wave <- 4

names(w5) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w5$wave <- 5

names(w6) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w6$wave <- 6

names(w7) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w7$wave <- 7

names(w8) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w8$wave <- 8

names(w9) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w9$wave <- 9

names(w10) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w10$wave <- 10

names(w11) <- c('id', 'X101', 'X103', 'X104', 'X105', 'X109', 'X110', 'X111', 'X112', 'X113', 'X114', 'X115', 'X116', 'X117', 'X118', 'X119', 'X120', 'X121', 'X122', 'X123', 'X124', 'X125', 'X126', 'X127', 'X128', 'X129', 'X130', 'X131', 'X132', 'X133', 'X134', 'X143', 'X144', 'X145', 'X146', 'X151', 'X152', 'X153', 'X154')

w11$wave <- 11

full <- rbind(w1, w2, w3, w4, w5, w6, w7, w8, w9, w10, w11)
#full <- full[,c(1, 44, 2:ncol(full))]
#gets rid of extra wave variable
#full <- full[,-45]

######
#selecting the background variable waves that correspond to my Politics and Values waves
dem$wave2 <- NA
dem$wave2 <- ifelse(dem$wave == 200712 | dem$wave == 200803, 1, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 200812 | dem$wave == 200901, 2, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 200912 | dem$wave == 201001 | dem$wave == 201002, 3, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201012 | dem$wave == 201101, 4, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201112 | dem$wave == 201201, 5, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201212 | dem$wave == 201301, 6, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201312 | dem$wave == 201401, 7, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201512 | dem$wave == 201601, 8, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201612 | dem$wave == 201701, 9, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201712 | dem$wave == 201801 | dem$wave == 201802 | dem$wave == 201803, 10, dem$wave2)
dem$wave2 <- ifelse(dem$wave == 201812 | dem$wave == 201901 | dem$wave == 201902 | dem$wave == 201903, 11, dem$wave2)

#getting rid of individuals younger than 18
dem$age <- NA
dem$age <- ifelse(dem$gebjaar <= 1990 & dem$wave2 == 1, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1991 & dem$wave2 == 2, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1992 & dem$wave2 == 3, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1993 & dem$wave2 == 4, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1994 & dem$wave2 == 5, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1995 & dem$wave2 == 6, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1996 & dem$wave2 == 7, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1997 & dem$wave2 == 8, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1998 & dem$wave2 == 9, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 1999 & dem$wave2 == 10, 1, dem$age)
dem$age <- ifelse(dem$gebjaar <= 2000 & dem$wave2 == 11, 1, dem$age)
#dem[is.na(dem)] <- 0
dem <- subset(dem, age == 1)
#turning birth year into age
dem$age <- NA
dem$age <- ifelse(dem$wave2 == 1, 2007-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 2, 2008-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 3, 2009-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 4, 2010-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 5, 2011-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 6, 2012-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 7, 2013-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 8, 2015-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 9, 2016-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 10, 2017-dem$gebjaar, dem$age)
dem$age <- ifelse(dem$wave2 == 11, 2018-dem$gebjaar, dem$age)

#cleaning marriage variable
dem$married <- dem$burgstat
dem$married <- ifelse(dem$burgstat == "Separated" | dem$burgstat == "Divorced" | dem$burgstat == "Widow or widower", "Previously married", dem$married)
dem$married <- ifelse(dem$married== 1, "Married", dem$married)
dem$married <- ifelse(dem$married== 5, "Never married", dem$married)

#cleaning urban variable
dem$urban <- dem$sted
dem$urban <- ifelse(dem$sted == "Extremely urban" | dem$sted == "Very urban", "Urban", dem$urban)
dem$urban <- ifelse(dem$urban == 3 | dem$urban == 4, "Suburban", dem$urban)
dem$urban <- ifelse(dem$urban == 5, "Rural", dem$urban)

#cleaning job variable
dem$job <- dem$belbezig
dem$job <- ifelse(dem$job == "Takes care of the housekeeping", "Homemaker", dem$job)
dem$job <- ifelse(dem$job == 1 | dem$job == 2 | dem$job == 3, "Employed", dem$job)
dem$job <- ifelse(dem$job == 7, "In school", dem$job)
dem$job <- ifelse(dem$job == 10 | dem$job == 12 | dem$job == 13 | dem$job == 14, "Other", dem$job)
dem$job <- ifelse(dem$job == 9, "Pensioner", dem$job)
dem$job <- ifelse(dem$job == 4 | dem$job == 5 | dem$job == 6 | dem$job == 11, "Unemployed", dem$job)

#cleaning income variable
dem$income <- as.numeric(dem$brutocat)
dem$income[dem$income == 14] <- 'dk'

#cleaning education variable
dem$education <- as.numeric(dem$oplmet)
dem$education <- ifelse(dem$education == 1 | dem$education == 2 | dem$education == 8 | dem$education == 9, "Less than HS", dem$education)
dem$education <- ifelse(dem$education == 3, "HS", dem$education)
dem$education <- ifelse(dem$education == 4 | dem$education == 5, "Vocational", dem$education)
dem$education <- ifelse(dem$education == 6, "College", dem$education)
dem$education <- ifelse(dem$education == 7, "Other", dem$education)

#getsrid of original education, etc variables that were recoded above
dem <- dem[,-c(5:9)]
#gets rid of background variable wave variable
dem <- dem[, -2]
#reorder columns
dem <- dem[,c(1, 4, 2, 3, 5:ncol(dem))]
names(dem)[c(1:2)] <- c('id', 'wave')

fullds <- merge(dem, full, by = c('id', 'wave'), all = F)

#individuals may be in the data more than once; this code removes duplicates
fullds$info_variable <- paste(fullds$geslacht, fullds$age, fullds$wave, fullds$married, fullds$urban, fullds$job, fullds$income, fullds$education, fullds$id, sep = ' ')
dups <- fullds[!duplicated(fullds$info_variable), ]

#getting rid of duplicated rows that are duplicated because they have NAs in the pasted variable
dups$eliminate <- ifelse(grepl('NA', dups$info_variable), 1, 0)
dups <- subset(dups, dups$eliminate == 0)

#finding which duplicated rows are remaining and removing them
dups$id2 <- paste(dups$id, dups$wave, sep = ' ')
dups$duplicated <- ifelse(duplicated(dups$id2), 1, 0)
dups <- subset(dups, duplicated == 0)
#removing info variables
dups <- dups[,-c(49:ncol(dups))]

#replacing the full dataset with the one that has removed duplicates
fullds <- dups

#Because there is a lot of missingness in the income variable, I am going to impute it later. The variable's NAs previously had a placeholder to avoid removing any rows with missing income in the previous steps. This line is now turning the placeholder into NA
fullds$income[fullds$income == 'dk'] <- NA

#####
names(abortion) <- c('id', 'abortion1', 'abortion2','abortion3', 'abortion4', 'abortion5', 'abortion6', 'abortion7', 'abortion8', 'abortion9', 'abortion10')

abortionlong <- reshape(abortion, varying = 2:11, sep = "", direction = 'long')

names(abortionlong)[2] <- 'wave'

full2 <- merge(fullds, abortionlong, by = c('id', 'wave'), all = F)

#######
#trying to remove rows with tons of NAs
full2 <- full2[!with(full2,is.na(X104) & is.na(X115)),]
full2 <- full2[!with(full2,is.na(X115) & is.na(X143)),]

full2$X101[full2$X101 == 999] <- NA
full2$X103[full2$X103 == 99] <- NA
full2$X104[full2$X104 == 99] <- NA
full2$X105[full2$X105 == 99] <- NA

full2$abortion <- as.character(full2$abortion)
full2$abortion <- ifelse(full2$abortion == 'yes', 1, full2$abortion)
full2$abortion <- ifelse(full2$abortion == 'maybe', 2, full2$abortion)
full2$abortion <- ifelse(full2$abortion == 'no', 3, full2$abortion)
full2$abortion <- ifelse(full2$abortion == "I don't know", NA, full2$abortion)

###merging in religion
names(religion) <- c('id', 'religion1', 'religion2', 'religion3', 'religion4', 'religion5', 'religion6', 'religion7', 'religion8', 'religion9', 'religion10')

religion <- reshape(religion, varying = 2:11, sep = "", direction = 'long')
names(religion)[2] <- 'wave'
religion <- na.omit(religion)

full2t <- merge(full2, religion, by = c('id', 'wave'), all = T)
full2t <- full2t[!with(full2t,is.na(geslacht) & is.na(gebjaar)),]
full2 <- full2t

#recoding so that higher values indicate frequent attendance
full2$religion <- car::recode(full2$religion, "1 = 7; 2 = 6; 3 = 5; 4 = 4; 5 = 3; 6 = 2; 7 = 1")

full2$religion[full2$religion == 99] <- NA

#reverse coding so that low values indicate liberalism
full2$X103 <- car::recode(full2$X103, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X109 <- car::recode(full2$X109, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X112 <- car::recode(full2$X112, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X114 <- car::recode(full2$X114, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X115 <- car::recode(full2$X115, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X116 <- car::recode(full2$X116, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X118 <- car::recode(full2$X118, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X119 <- car::recode(full2$X119, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X122 <- car::recode(full2$X122, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X126 <- car::recode(full2$X126, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X127 <- car::recode(full2$X127, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X128 <- car::recode(full2$X128, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X129 <- car::recode(full2$X129, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
full2$X130 <- car::recode(full2$X130, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")

full2$income <- as.numeric(full2$income)
full2$abortion <- as.numeric(full2$abortion)

##### now I am imputing missing values
library(mice)
ini <- mice(full2, maxit = 0)
#this shows the imputation method for each variable. All should be "pmm". Variables with no method have no missingness and thus will not be imputed
meth <- ini$meth
meth 
pred <- ini$predictorMatrix 

#this code will take several hours to run
imp1 <- mice(full2, meth = meth, pred = pred, print = FALSE, nnet.MaxNWts = 5000)

#you may wish to save the imputation: saveRDS(imp1, 'imp1.RDS')

#this turns the imputation into a useable dataset
dat <- mice::complete(imp1)

#checking cronbach's alpha for gender
psych::alpha(dat[c(16:22, 31:37, 42:49, 54)])

library(mirt)
#making gender attitudes variable
mirt1 <- mirt(data=dat[c(16:22, 31:37, 42:49, 54)], model = 1, itemtype="graded", SE=TRUE, verbose=T, na.rm = T)

#extracting attitude for each individual
dat$fullgen <- fscores(mirt1)
dat$fullgen <- as.numeric(dat$fullgen)

#immigration attitude
psych::alpha(dat[c(14, 23:30)], check.keys = T)
#immigration alpha indicates two variables are reverse coded:
dat$X117<- car::recode(dat$X117, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")
dat$X121<- car::recode(dat$X121, "1 = 5; 2 = 4; 3 = 3; 4 = 2; 5 = 1")

mirt2 <- mirt(data=dat[c(14, 23:30)], model = 1, itemtype="graded", SE=TRUE, verbose=T, na.rm = T)

dat$im <- fscores(mirt2)
dat$im <- as.numeric(dat$im)

#making cohort variable
dat$cohort <- NA
dat$cohort <- ifelse(dat$gebjaar <=1949, 1, dat$cohort)
dat$cohort <- ifelse(dat$gebjaar >=1950 & dat$gebjaar <=1959, 2, dat$cohort)
dat$cohort <- ifelse(dat$gebjaar >=1960 & dat$gebjaar <=1969, 3, dat$cohort)
dat$cohort <- ifelse(dat$gebjaar >=1970 & dat$gebjaar <=1979, 4, dat$cohort)
dat$cohort <- ifelse(dat$gebjaar >=1980, 5, dat$cohort)
dat$cohort <- as.factor(dat$cohort)

#making cohort dummy
dat$cohort2 <- NA
dat$cohort2 <- ifelse(dat$gebjaar <=1969, 0, 1)
dat$cohort2 <- as.factor(dat$cohort2)

#save data to import into stata for further analyses
write.csv(dat, 'dat.csv')

###
#Figure 1
library(dplyr)
plotdat <- data.frame(dat$fullgen, dat$im, dat$cohort)
colnames(plotdat) <- c('fullgen', 'im', 'cohort')

#this creates a dummy variable that takes on a 1 for those with gender attitudes that are less conservative than average and immigration attitudes that are more opposed than average (liberal nationalists)
plotdat$freq <- ifelse(plotdat$fullgen < 0 & plotdat$im > 0, 1, 0)

#counting how many are in each cohort
plyr::count(plotdat$cohort)
#counting how many liberal nationalists in each cohort
plotdat %>% group_by(cohort) %>% count(freq)

#these are the calculated prevalences by cohort
#cohort 1: 2133/16510 = 0.1291944
#cohort 2: 2159/12921 = 0.1670923
#cohort 3: 2378/11976 = 0.1985638
#cohort 4: 2216/9669 = 0.2291861
#cohort 5: 2318/11454 = 0.2023747

#this creates a dummy variable that takes on a 1 for those with gender attitudes that are less conservative than average and immigration attitudes that are less opposed than average (liberals)
plotdat$freq2 <- ifelse(plotdat$fullgen < 0 & plotdat$im < 0 , 1, 0)

#counting how many liberals in each cohort
plotdat %>% group_by(cohort) %>% count(freq2)

#these are the calculated prevalences by cohort
#3768/16510 = 0.2282253
#4220/12921 = 0.3266001
#3717/11976 = 0.3103707
#2930/9669 = 0.3030303
#3797/11454 = 0.3314999

#this creates a dummy variable that takes on a 1 for those with gender attitudes that are more conservative than average and immigration attitudes that are more opposed than average (conservatives)
plotdat$freq3 <- ifelse(plotdat$fullgen > 0 & plotdat$im > 0 , 1, 0)

#counting how many conservatives in each cohort
plotdat %>% group_by(cohort) %>% count(freq3)

#these are the calculated prevalences by cohort
#6444/16510 = 0.3903089
#3812/12921 = 0.2950236
#3459/11976 = 0.2888277
#2485/9669 = 0.2570069
#2608/11454 = 0.2276934

#combining the prevalences into a dataframe
plotdat <- data.frame(cohort = 1:5, libnat = c(0.1291944, 0.1670923, 0.1985638, 0.2291861, 0.2023747), libs = c(0.2282253, 0.3266001, 0.3103707, 0.3030303, 0.3314999), cons = c(0.3903089, 0.2950236, 0.2888277, 0.2570069, 0.2276934))

plotdat <- plotdat %>%
  tidyr::gather(`libnat`:`cons`, key="Type", value="Mean")

library(ggplot2)
fig1 <- ggplot(plotdat, aes(y=Mean, x=cohort, color = Type)) + 
  geom_line() + geom_point() + scale_color_manual(name="Attitude Type", labels=c("Conservative",'Liberal Nationalist',"Liberal"), values = c("grey20", "grey50", 'grey70')) + labs(x = "Cohort", y = "Percentage") +theme_bw() + xlim('<1950', '1950-9', '1960-9', '1970-9', '1980+') + ylim(.1, .5)

#Figure 7
####
#first, finding respondents in wave 1 who were born between certain years
datsub1 <- subset(dat, wave == 1)

#1980-1985
datsub1$try <- ifelse(datsub1$gebjaar >= 1980 & datsub1$gebjaar <=1985, 1, NA)
#1970-1975
datsub1$try <- ifelse(datsub1$gebjaar >= 1970 & datsub1$gebjaar <=1975, 2, datsub1$try)
#1950-1955
datsub1$try <- ifelse(datsub1$gebjaar >= 1950 & datsub1$gebjaar <=1955, 3, datsub1$try)

#removing all other respondents
datsub1 <- subset(datsub1, !is.na(try))

#now finding which of these respondents are also present in wave 11
datsub11 <- subset(dat, wave == 11)
datsub11$try <- NA
datbind <- rbind(datsub1, datsub11)
res <- datbind %>% count(id)
res <- subset(res, n == 2)
dats <- merge(datbind, res, by = 'id')

#now filling in "try" for wave 11
dats$try <- ifelse(dats$wave == 11 & dats$gebjaar >= 1980 & dats$gebjaar <=1985, 1, dats$try)
dats$try <- ifelse(dats$wave == 11 & dats$gebjaar >= 1970 & dats$gebjaar <=1975, 2, dats$try)
dats$try <- ifelse(dats$wave == 11 & dats$gebjaar >= 1950 & dats$gebjaar <=1955, 3, dats$try)

#creating separate datasets for each group in each wave
dats1.1 <- subset(dats, try == 1 & wave == 1)
dats1.11 <- subset(dats, try == 1 & wave == 11)
dats2.1 <- subset(dats, try == 2 & wave == 1)
dats2.11 <- subset(dats, try == 2 & wave == 11)
dats3.1 <- subset(dats, try == 3 & wave == 1)
dats3.11 <- subset(dats, try == 3 & wave == 11)

#running correlation tests on each subset and extracting correlation and confidence intervals
cor1.1 <- cor.test(dats1.1$fullgen, dats1.1$im)$estimate
ci1.1 <- cor.test(dats1.1$fullgen, dats1.1$im)$conf.int
cor1.11 <- cor.test(dats1.11$fullgen, dats1.11$im)$estimate
ci1.11 <- cor.test(dats1.11$fullgen, dats1.11$im)$conf.int

cor2.1 <- cor.test(dats2.1$fullgen, dats2.1$im)$estimate
ci2.1 <- cor.test(dats2.1$fullgen, dats2.1$im)$conf.int
cor2.11 <- cor.test(dats2.11$fullgen, dats2.11$im)$estimate
ci2.11 <- cor.test(dats2.11$fullgen, dats2.11$im)$conf.int

cor3.1 <- cor.test(dats3.1$fullgen, dats3.1$im)$estimate
ci3.1 <- cor.test(dats3.1$fullgen, dats3.1$im)$conf.int
cor3.11 <- cor.test(dats3.11$fullgen, dats3.11$im)$estimate
ci3.11 <- cor.test(dats3.11$fullgen, dats3.11$im)$conf.int

try <- rbind(cor1.1, cor1.11, cor2.1, cor2.11, cor3.1, cor3.11)
try2 <- rbind(ci1.1[1], ci1.11[1], ci2.1[1], ci2.11[1], ci3.1[1], ci3.11[1])
try3 <- rbind(ci1.1[2], ci1.11[2], ci2.1[2], ci2.11[2], ci3.1[2], ci3.11[2])

#binding into a dataframe
cors <- as.data.frame(cbind(try, try2, try3))
cors$wave <- c('2007-8', '2018-9', '2007-8', '2018-9', '2007-8', '2018-9')
cors$wave <- as.factor(cors$wave)
cors$cohort <- c('1980-85', '1980-85', '1970-75', '1970-75', '1950-55', '1950-55')
cors$cohort <- as.factor(cors$cohort)

fig7 <- ggplot() + geom_point(aes(cors$wave, cors$cor, color = cors$cohort), size = 2.25, position=position_dodge(width=0.5))+ scale_color_manual(values = c('1980-85' = 'grey5', '1970-75' = 'grey45', '1950-55' = 'grey70'))+ labs(y = 'Gender and Immigration Attitudes Correlation', x = 'Wave', color = 'Birth Year') + geom_linerange(aes(x = cors$wave, ymin = cors$V2, ymax = cors$V3, color = cors$cohort), position=position_dodge(width=0.5)) + scale_color_manual(values = c('1980-85' = 'grey5', '1970-75' = 'grey45', '1950-55' = 'grey70')) + geom_hline(aes(yintercept = 0), linetype = 'dashed') + theme_bw()

